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PREFACE 


Bounded  nominal  paths  can  be  constructed  in  the  vicinity  of  the 
interior  equilibrium  point  (sometimes  called  a  libration  or  Lagrange 
point)  for  the  Sun-Ear th+Moon  Elliptic  Restricted  Three-Body  Problem. 
The  nominal  paths  may  be  computed  using  an  analytic  approximation  or, 
more  accurately,  using  a  numerical  integration  routine.  Numerical 
integration  is  used  to  generate  the  periodic  or  quasi-periodic 
reference  trajectories  in  this  effort.  The  output  of  the  routine  will 
be  numerical  values  for  each  of  the  six  states  (three  position  and 
three  velocity)  at  each  of  the  integration  time  steps.  Linearization 
of  both  the  equations  of  motion  and  of  the  equations  representing  the 
tracking  solution  assumes  access  to  a  continuous  representation  of  the 
spacecraft’s  orbit.  Follow-on  research  that  investigates  tracking 
errors  or  station-keeping  costs  may  also  need  a  continuous  (and  smooth) 
representation  of  the  six  states  or,  at  least,  may  need  access  to  an 
interpolation  routine.  Consequently,  this  work  explores  the  generation 
of  curves  through  the  numerical  data  representing  the  libration  point 
orbits  in  the  Sun-Ear th+Moon  ER3BP;  these  orbits  are  computed  in  the 
vicinity  of  the  interior  libration  point  between  the  Sun  and  the 
Earth/Moon  barycenter.  This  effort  is  supported  by  the  Frank  J.  Seiler 
Research  Laboratory  and  has  been  conducted  as  doctoral  research  under 
the  direction  of  Professor  K.C.  Howell,  School  of  Aeronautics  and 
Astronautics,  Purdue  University,  West  Lafayette,  Indiana. 
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INTRODUCTION 


With  the  expansion  of  space  exploration  programs  worldwide, 
interest  has  increased  in  the  design  of  innovative,  complex,  and  yet 
low-cost  spacecraft  trajectories  that  meet  demanding  mission 
requirements.  In  most  of  the  missions  flown  in  the  last  few  decades, 
the  spacecraft  spent  the  majority  of  the  flight  time  in  a  force 
environment  dominated  by  a  single  gravitational  field.  For  the 
preliminary  mission  analysis  in  these  cases,  additional  attracting 
bodies  and  other  forces  could  be  modeled,  when  required,  as  perturbing 
influences.  Analysis  of  some  recently  proposed  and  more  adventurous 
missions,  such  as  those  involving  libration  point  orbits,  will  require 
dynamic  models  of  higher  complexity,  since  at  least  two  gravitational 
fields  are  of  nearly  equal  influence  on  the  spacecraft  throughout  the 
majority  of  the  mission.  Thus,  trajectories  determined  for  a  system 
consisting  of  numerous  gravitational  forces  have  been  of  particular 
theoretical  and  practical  interest  in  recent  years. 

One  type  of  many-body  problem,  motion  within  a  three-body  system 
of  forces,  has  a  wide  range  of  applications.  The  general  problem  of 
three  bodies  assumes  that  each  body  has  finite  mass  and  that  the  motion 
is  a  result  of  mutual  gravitational  attraction.  When  the  mass  of  one 
of  the  three  bodies  is  assumed  to  be  sufficiently  small  (infinitesimal) 
so  that  it  does  not  affect  the  motion  of  the  other  two  bodies 
(primaries)  in  the  system,  the  "restricted  three-body  problem"  results. 
The  primaries  can  be  further  assumed  to  be  moving  in  known  elliptic  or 
circular  orbits  about  their  common  center  of  mass.  Therefore,  the 
elliptic  restricted  three-body  problem,  where  the  primaries  are  assumed 
to  be  in  known  elliptic  orbits,  may  be  considered  a  reasonably 
approximate  model  for  a  spacecraft  moving  within  the  gravitational 
fields  of  the  Sun  and  the  Earth,  for  instance. 
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In  the  formulation  of  the  restricted  three-body  problem,  one  mass 
is  defined  as  infinitesimal  relative  to  the  remaining  two  masses 
(primaries).  The  primaries,  unaffected  by  the  infinitesimal  mass,  move 
under  their  mutual  gravitational  attractions.  In  the  elliptic 
restricted  three-body  problem  (ER3BP),  the  primaries  are  assumed  to 
move  on  elliptic  paths.  If  the  eccentricity  of  the  primaries1  orbit  is 
assumed  to  be  zero,  the  circular  restricted  three-body  problem  (CR3BP) 
results.  Even  for  known  primary  motion,  however,  a  general, 
closed-form  solution  for  motion  of  the  third  body  of  infinitesimal  mass 
does  not  exist.  In  the  restricted  three-body  problem  (ER3BP  or  CR3BP), 
five  equilibrium  (libration)  solutions  can  be  found.  These  equilibrium 
points,  sometimes  also  called  Lagrange  points,  are  particular  solutions 
of  the  equations  of  motion  governing  the  path  of  the  infinitesimal  mass 
moving  within  the  gravitational  fields  of  the  primaries. 

The  equilibrium  points  are  defined  relative  to  a  coordinate  system 
rotating  with  the  primaries.  At  these  locations,  the  forces  on  the 
spacecraft  are  in  equilibrium.  These  forces  include  the  gravitational 
forces  from  the  massive  bodies  and  the  centrifugal  force  associated 
with  the  rotation  of  the  system.  (The  addition  of  solar  radiation 
pressure  to  the  force  model  changes  the  locations  of  the  five  Lagrange 
points,  although  they  can  still  be  defined,  and  these  solar  radiation 
effects  are  discussed  in  Gordon115.)  The  libration  points  are  located 
in  the  plane  of  primary  rotation.  Three  of  the  libration  points  are  on 
the  line  between  the  two  massive  bodies,  and  one  of  these  collinear 
points  is  interior  to  the  primaries.  The  last  two  points  are  at  the 
vertices  of  two  equilateral  triangles  in  the  plane  of  primary  rotation. 
The  triangles  have  a  common  base  that  is  the  line  between  the  primary 
masses. 

For  the  CR3BP,  the  five  libration  points  are  stationary  relative 
to  the  rotating  reference  frame.  If  the  problem  is  generalized  to  the 
ER3BP,  the  libration  points  pulsate  as  the  distance  between  the 
primaries  varies  with  time.  In  both  the  circular  and  elliptic 
restricted  problems,  two-dimensional  and  three-dimensional 
trajectories,  both  periodic  and  quasi-periodic  paths,  can  be  computed 
in  the  vicinity  of  these  libration  points. 
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Three-dimensional,  periodic  "halo"  orbits  in  the  vicinity  of  the 

collinear  libration  points  have  been  studied  since  the  late  1960s. 

Early  work  concerning  these  orbits  was  motivated  by  studies  related  to 

exploring  the  far  side  of  the  Moon.  These  studies  were  completed  in 

support  of  the  planned  Apollo  18  lunar  exploration  mission  that  was 

later  canceled.  Robert  Farquhar  coined  the  term  "halo"  to  describe  a 

three-dimensional,  periodic  orbit  near  a  libration  point  on  the  far 

(21 

side  of  the  Moon  in  the  Earth-Moon  system.  These  orbits  were 

designed  to  be  large  enough  so  that  the  spacecraft  would  be  constantly 
in  view  of  the  Earth  and  would  thus  appear  as  a  halo  around  the  Moon. 
A  communications  station  in  this  type  of  orbit  could  maintain  constant 
contact  between  the  Earth  and  a  lunar  experimentation  station  on  the 
far  side  of  the  Moon. *31 

Quasi-periodic  orbits  near  libration  points  are  also  currently  of 

great  research  interest.  The  variations  in  size  and  shape  that  a 

quasi-periodic  orbit  can  exhibit  may  add  valuable  flexibility  for 

mission  planning.  This  type  of  bounded,  three-dimensional  libration 

point  trajectory  is  called  a  L-issajous  orbit  since  specific  planar 

projections  of  these  quasi-periodic  trajectories  may  look  like  a 

special  type  of  "Lissajous"  curve.  Physicist  Jules  Antoine  Lissajous 

(1822-1880)  investigated  curves  that  were  generated  by  compounding 

simple  harmonic  motions  at  right  angles,  and  he  delivered  a  paper  >n 

this  subject  to  the  Paris  Academy  of  Sciences  in  1857.  Nathaniel 

Bowditch  of  Salem,  Massachusetts,  had  conducted  some  similar  work  in 

1815.  Lissajous  curves  have  a  wide  variety  of  shapes  that  depend  on 

the  frequency,  phase,  and  amplitude  of  the  orthogonal  components  of  the 
[4  S] 

motion.  ’  When  the  in-plane  and  the  (orthogonal)  out-of-plane 
frequencies  of  the  linearized  motion  are  nearly  (but  not)  equal,  the 
resulting  path  is  typically  called  a  Lissajous  trajectory. 

A  method  to  generate  approximations  for  this  type  of 
quasi-periodic  orbital  path  was  developed  analytically  by  Farquhar  and 
Kamel  in  1973. 161  They  derived  a  third-order  approximate  analytic 
solution  for  a  translunar  libration  point  orbit  in  the  Earth-Moon  ER3BP 
that  also  included  solar  gravity  perturbations.  In  1975,  Richardson 
and  Cary  then  developed  a  fourth-order  analytic  Lissajous  approximation 
in  the  Sun-Earth+Moon  barycenter  system.  The  notation  "Earth+Moon" 


3 


indicates  that  the  Earth  and  the  Moon  are  treated  as  one  body  with  mass 

center  at  the  Earth-Moon  .arycenter.  In  consideration  of  the  lunar 

perturbation,  Farquhar  has  shown  that  the  accuracy  of  solutions  in  the 

Sun-Earth  CR3BP  can  be  enhanced  if  the  collinear  libration  points  are 

defined  along  the  line  between  the  Sun  and  the  center  of  mass  of  the 

Earth  and  the  Moon.  Since  1975,  a  few  researchers  have  considered 

methods  to  numerically  generate  Llssajous  trajectories,  but  the  lack  of 

periodicity  of  a  Llssajous  path  complicates  numerical  construction  of 

bounded  trajectories.  Howell  and  Pernicka  have  developed  a  numerical 

technique  for  determination  of  three-dimensional,  bounded  Llssajous 

[9-14] 

trajectories  of  arbitrary  size  and  duration. 

The  numerical  integration  routine  computes  bounded  periodic  or 
quasi-periodic  libration  point  orbits.  The  numerical  ouput  (three 
position  and  three  velocity  states)  is  indexed  by  integration  time 
steps;  follow-on  research  may  require  a  continuous  representation 
of  the  orbit  instead  of  a  tabular  listing  of  the  numerical  data.  The 
first  chapter  briefly  covers  the  background  of  the  p-oblem.  The 
following  chapter  then  discusses  the  four  methods  of  curve  fitting 
investigated  and  the  method  selected  to  be  used  in  follow-on  research. 
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CHAPTER  1:  BACKGROUND 


In  this  chapter,  the  elliptic  restricted  three-body  problem  and 
the  associated  coordinate  systems  are  reviewed;  the  equations  of  motion 
for  an  infinitesimal  mass  moving  in  the  gravity  fields  of  two  massive 
bodies  are  then  presented.  Next,  locations  of  the  libration  points  are 
discussed.  Finally,  the  state  transition  matrix  and  the  construction 
of  bounded  nominal  orbits  near  the  collinear  Lagrange  points  are 
summarized.  A  more  thorough  discussion  of  these  topics  is  presented  in 
Gordon. 111 


A.  Elliptic  Restricted  Three-Body  Problem 

The  elliptic  restricted  three-body  problem  is  a  simplification  of 
the  general  problem  of  three  bodies.  In  the  general  three-body 
problem,  each  of  the  three  bodies  is  assumed  to  be  a  particle  of  finite 
mass  and,  thus,  exerts  an  influence  on  the  motion  of  each  of  the  other 
bodies.  Neither  the  general  nor  the  restricted  problem  of  three  bodies 
has  a  general  closed-form  solution.  However,  when  problem 
simplifications  are  made,  particular  solutions  can  be  determined.  If 
the  mass  of  one  of  the  bodies  is  restricted  to  be  infinitesimal,  such 
that  it  does  not  affect  the  motion  of  the  other  two  massive  bodies 
(primaries),  the  restricted  three-body  model  results.  The  primaries 
are  assumed  to  be  in  known  elliptic  (ER3BP)  or  circular  (CR3BP)  orbits 
about  their  common  mass  center  (barycenter) .  The  problem  can  then  be 
completely  described  by  a  single  second-order  vector  differential 
equation  with  variables  appropriately  defined  for  a  specified 
coordinate  frame. 
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B.  Coordinate  Systems 


The  two  standard  coordinate  systems  used  in  the  analysis  of  this 

problem  have  a  common  origin  at  the  center  of  mass  (barycenter)  of  the 

primaries.  Primaries  with  masses  m  and  m  such  that  m  s  m  are 

1212 

assumed  here,  although  this  distinction  is  arbitrary.  The 

infinitesimal  mass  is  denoted  as  m  .  These  masses  (m  ,m  ,m  ) 

3  12  3 

correspond  to  particles  situated  at  points  P  ,  P  ,  and  P  , 

12  3 

respectively.  The  barycenter  is  denoted  as  "B,"  and  the  resulting 

arrangement  is  shown  in  Figure  1-1.  The  rotating  coordinate  system  is 

defined  as  x  y  z  ,  and  the  inertial  system  is  identified  as  x  y  z  . 

ft  R  R  ill 

Note  that  both  coordinate  systems  are  right-handed,  and  the  x  and  y 
axes  for  both  systems  are  in  the  plane  of  motion  of  the  primaries.  The 
Xj  axis  is,  of  course,  assumed  to  be  oriented  in  some  fixed  direction; 
in  this  specific  formulation  of  the  problem,  it  is  assumed  to  be 
parallel  to  a  vector  defined  with  a  base  point  at  the  Sun  and  directed 
toward  periapsis  of  the  Earth’s  orbit.  The  rotating  x  axis  is  defined 

R 

along  the  line  that  joins  the  primaries  and  is  directed  from  the  larger 

toward  the  smaller  primary.  The  z  axes  are  coincident  and  are  directed 

parallel  to  the  primary  system  angular  momentum  vector.  The  y  axis 

R 

completes  the  right-handed  x  y  z  system. 

R  R  R 


C.  Equations  of  Motion 

Newtonian  mechanics  are  used  to  formulate  the  equations  of  motion 
for  m3  (the  spacecraft)  relative  to  B  as  observed  in  the  inertial 
reference  frame.  The  sum  of  the  forces  on  m3  resulting  from  both  the 
gravity  fields  of  masses  m^  (the  Sun)  and  m2  (the  Earth-Moon 
barycenter)  and  from  the  solar  radiation  pressure  can  be  used  to 
produce  the  following  second-order  vector  differential  equation: 


mm* 

P 


m 

=  -  G  (-i-)  d  -  G  ( 
d3 


)  r  +  ( 


)  d. 


(1-1) 
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The  overbar  denotes  a  vector,  and  primes  indicate  differentiation 
with  respect  to  dimensional  time.  All  quantities  are  dimensional,  as 
appropriate,  and  the  quantity  "G"  is  the  universal  gravitiational 
constant.  The  scalars  "d"  and  11  r"  lr  equation  (1-1)  denote  the 
magnitudes  of  the  vectors  d  and  r,  respectively,  as  depicted  in 
Figure  1-1.  The  dimensionless  scalar  “k"  is  the  solar  reflectivity 
constant,  and  "S"  is  the  solar  radiation  pressure  constant.  The 
formulation  of  the  solar  radiation  force  model  and  the  values  for  the 
solar  radiation  constants  are  derived  from  previous  work  by  Bell.1151 
The  numerical  values  for  these  constants  are  selected  from 
characteristic  data  for  a  spacecraft  previously  in  a  libration  point 
orbit: 


S  = 


A  S  r 
o  o 

c  m  ’ 


(1-2) 


where  k  =  1.2561,  S=  1352.098  kg/sec3,  r  =  1.4959787x10s  km,  c  =  the 

speed  of  light  =  2.998x10  m/sec,  A  =  surface  area  of  the  spacecraft’s 

sun-facing  side  =  3.55  m2,  and  m3  =  mass  of  the  spacecraft  =  435  kg. 

The  constant  Sq  is  the  solar  light  flux  measured  at  one  astronomical 

unit  (A.U.  ) ,  and  rQ  is  the  nominal  distance  associated  with  the 

measurement  of  the  solar  flux  S  .  [15]  The  value  of  A,  also  termed  the 

o 

effective  cross-sectional  area  of  the  spacecraft,  is  considered  to  be 
constant  in  this  work.  The  constant  k  is  a  material  parameter 
dependent  on  the  absorptivity  of  the  spacecraft  surface  and  is 
generally  confined  to  a  range  of  0  i  k  s  2.0.I1S) 

The  position  vector  p  is  written  in  rotating  components  as 


p=xx+  yy+zz  (1-3) 

K  r  1  7r  r 


where  x  ,y  ,z  are  unit  vectors.  The  velocity  and  the  acceleration  of 

R  R  R 

the  spacecraft  (particle  with  mass  m3)  relative  to  the  barycenter  B 
as  observed  in  the  inertial  reference  frame  can  then  be  derived.  The 
following  kinematic  expression  for  p"  can  be  derived: 
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(1-4) 


p"  =  (x"-e"y-2e,y,-e,2x)x  +(y"+0"x+20'x' -0,2y)y  +z"z  . 

H  R  R 

Notice  that  equation  (1-4)  includes  terms  with  the  angular 
velocity,  0'(t),  and  angular  acceleration,  0"(t),  of  the  rotating 
reference  frame  with  respect  to  the  inertial  frame.  By  using  known 
elliptic  orbital  elements,  readily  computed  expressions  for  0'  and  0" 
can  be  found: 


Vatl-e) (l+e)  /G(m  +m  ) 
_ 1  2 

a2  (1  -  e  cos(E))2 


„  -2G(m  +m  )/(”l-e)(l+e)  e  sin(E) 

0  =  - - — - 

a3(l  -  e  cos(E))4 


(1-6) 


where  e  is  the  eccentricity,  E  is  the  eccentric  anomaly,  and  a  is 
the  semi-major  axis  of  the  primary  orbit. 

Three  scalar  equations  of  motion  for  P3  can  be  derived  using  the 
dimensional  equations  (1-1),  (1-2),  (1-4),  (1-5),  and  (1-6);  however, 
for  convenience,  the  following  scaling  factors  are  typically 
introduced: 


(1)  The  sum  of  the  masses  of  the  primaries  equals  one 

mass  unit,  (m  +  m  =1  unit  of  mass) 

1  2 

(2)  The  mean  distance  between  the  primaries  equals  one 
unit  of  distance,  (a  =  1  unit  of  distance) 

(3)  The  universal  gravitational  constant  is  equal  to  one 

unit  by  proper  choice  of  characteristic  time, 
(characteristic  time  =  1/n  ;  thus  G  =  1  unit) 

mean 


The  dimensional  equations  of  motion  can  be  simplified  and  scaled 
by  introducing  the  characteristic  quantities  defined  above  and  by 
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introducing  the  nondimensional  mass  ratio  p,  "psuedo-potential"  U,  and 
the  scaled  solar  radiation  constant  s: 


and 


M  = 


m  +  m 
1  2 


(1-7) 


U  = 


(1-p) 

d 


+ 


•2,  2 

0  (x 


+ 


2 

y 


k  s 
d 


(1-8) 


where  the  dot  denotes  the  derivative  with  respect  to  characteristic 
time.  The  scaled  solar  radiation  constant,  s,  is  derived  from  the 
dimensional  radiation  constant  denoted  as  S  in  equation  (1-2)  and,  by 
using  the  characteristic  quantities  described  above,  its  value  is 
calculated  as  s  =  6. 206597029461384xl0-6  nondimensional  units.:  Then, 

the  vector  magnitudes,  "d"  and  "r,"  are  written  in  terms  of  scaled 
quantities  as: 


,  r  r  rx  \2  2  2,1/2 

d  =  [  (x  +  |i  R)  +  y  +  z) 


(1-9) 


r  f  rt.  n\2  2  2  ,  1/2 

r  =  [ (x  -  R  +  p  R)  +y  +  z  ] 


(1-10) 


The  three  scalar  second-order  differential  equations  that  result 
can  be  written  in  terms  of  characteristic  quantities  as 
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au 

ax 


(i-ii) 


x  -  2  0  y  = 


+  0  y  =  U  +  0  y, 


y  +  2  0  x 


au 

ay 


0  x  =  U  -  0  x, 
y 


(1-12) 


(1-13' 


The  scaled  equations  for  the  angular  velocity  and  angular 
acceleration  of  the  rotating  frame  relative  to  the  inertial  frame  also 
simplify: 


V(l-e) (1+e) 

9  =  - ,  (1-14) 

(1  -  e  cos(E) ) 


-2y/[\-e)  (1+e)  e  sin(E) 

0  =  - - - .  (1-15) 

(1  -  e  cos(E)j 


If  the  primaries  are  assumed  to  be  moving  in  a  circular  orbit, 
then  0  =  0,  R  =  0  =  1,  and  equations  (1-11),  (1-12),  (1-13)  reduce  to 
three  scalar  equations  in  the  simplified  form: 
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(1-16) 


x  -  2  y 


y  +  2 


(1-17) 


au 

az 


(1-18) 


The  scalar  equations  (1-11),  (1-12),  and  (1-13)  corresponding  to 
the  elliptic  restricted  problem  or  equations  (1-16),  (1-17),  and  (1-18) 
derived  for  the  circular  restricted  problem  can  be  used  to  locate  the 
five  libration  points  in  the  rotating  reference  frame. 


D.  Locations  of  the  Lagrangian  (Libration)  Points 

By  using  scalar  equations  (1-16),  (1-17),  and  (1-18)  for  motion  in 
the  CR3BP,  the  locations  of  the  stationary  equilibrium  points  can  be 
determined.  Equations  (1-11),  (1-12),  and  (1-13)  can  be  used  to 
determine  ratios  of  distances  that  are  constant  in  the  ER3BP;  these 
ratios  are  related  to  the  locations  of  libration  points  that  have  been 
defined  in  the  ER3BP  and  that  "pulsate"  with  respect  to  the  rotating 
reference  frame  as  the  distance  between  the  primaries  varies  with  time. 


1 .  The  CR3BP 

In  the  CR3BP,  the  five  libration  points  are  equilibrium  points  and 
arc  stationary  with  respect  to  the  rotating  coordinate  frame,  that  is, 
they  are  locations  at  which  the  foj  ces  on  the  third  body  sum  to  zero. 
The  arrangement  of  points  and  the  corresponding  nondimensional 
distances  are  depicted  in  Figure  1-2.  Note  that  three  of  the  libration 
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points  (L  ,  L2>  L3)  are  collinear  with  the  primaries;  one  collinear 
point  ( )  is  interior  to  the  primaries.  The  remaining  two  points  (L^ 
and  Lg)  are  located  at  the  vertices  of  two  equilateral  triangles  that 
are  in  the  plane  of  primary  rotation  and  that  have  a  common  base 
between  the  primaries. 

In  the  CR3BP,  the  libration  points  are  stationary  in  the  rotating 
coordinate  frame.  Stationary  points  are  defined  as  points  at  which  the 
relative  velocity  and  acceleration  are  zero,  such  that 


x=y=z=x=y=z=0.  (1-19) 

By  using  equations  (1-19)  in  equations  (1-16)  through  (1-18),  the 
useful  conditions  U  =  U  =  U  =  0  are  found.  The  three  collinear 

x  y  z 

libration  points  can  be  readily  located  by  further  noting  that  y  =  z  = 

0  for  the  points  located  on  the  rotating  x  axis. 

R 


2.  The  ER3BP 

Five  libration  points  also  exist  in  the  ER3BP,  but  they  are  not 

stationary  relative  to  the  rotating  frame;  rather,  the  collinear  points 

pulsate  along  the  x  axis,  and  the  triangular  points  pulsate  relative 

R 

to  both  the  x  and  the  y  axes  as  the  distance  between  the  primaries 

R  R 

varies  v/ith  time.  The  equilibrium  solutions  can  be  located  by  using 
equations  (1-11)  through  (1-13)  to  find  ratios  of  certain  distances 
that  are,  in  fact,  constant  in  the  problem.  The  collinear  libration 
points  in  the  ER3BP  can  be  found  by  assuming  x  *  0,  x  *  0,  and  y  =  y  = 
z=z=y=z=0.  The  relative  locations  of  the  libration  points  in 
the  ER3BP  are  depicted  in  Figure  1-3. 
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E.  State  Transition  Matrix 


The  state  transition  matrix  is  used  in  the  calculation  of  the 
acceptable  nominal  trajectory,  and  it  must  also  be  available  at  varying 
time  intervals  along  the  nominal  path  for  orbit  determination  error 
analysis  investigations  and  station-keeping  studies.  The  transition 
matrix  is  derived  in  connection  with  a  linearizing  analysis. 

The  equations  of  motion  for  the  infinitesimal  mass  in  the  ER3BP 
can  be  linearized  about  a  reference  trajectory  (nominal  path)  that  is  a 
solution  of  the  differential  equations.  The  states,  three  position  and 
three  velocity,  and  the  state  vector  x  are  defined  as 


x 


l 


x,  x  =  y,  x  =  z 

2  7  3 


X 


6 


Z, 


(1-20) 


and 


x  =  (x. 


*6] 


(1-21) 


The  reference  trajectory  is  defined  as  x  .  Therefore,  using  a 

REF 

Taylor’s  series  approach,  the  expansion  about  the  reference  path  is 
written  in  the  form  of  the  first-order  variational  equation 


— (  x)  =  x  =  A ( t )  x 
dt 


(1-22) 


where  x  =  x  -  x  is  understood  to  be  the  vector  of  residuals  relative 

REF 

to  the  nominal  solution,  and  the  matrix  A ( t )  contains  the  first-order 
terms  in  the  Taylor’s  series  expansion  of  the  equations  of  motion  about 
the  nominal  or  reference  solution  of  interest. 
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Using  equations  (1-11)  through  (1-13),  A(t)  can  be  expressed  as 


A(t)  = 


0 

Urr+80 


I 

2en 


(1-23) 


where  all  four  submatrices  are  dimension  3x3  and 


with 


Uxx 

Uxy 

Uxz 

Uyx 

Uyy 

Uyz 

Uzx 

Uzy 

Uzz 

1 

0 

0 


O' 

0 

0 


(1-24) 


In  equation  (1-24),  the  notation  is  simplified  for  the  partial 
derivatives;  for  instance 


S2U 


dx 


2 


Uxx. 


The  matrix  A(t)  can  then  be  evaluated  along  the  reference  trajectory. 

The  vector  differential  equation  (1-22)  governing  the  state 
variations  from  the  nominal  path  has  a  solution  of  the  form 
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(1-25) 


x(t)  =  *(t,t  )  x(t  ) 

0  0 

where  *(t,tQ)  is  the  state  transition  matrix  at  time  "t"  relative  to 
time  "t  . "  The  matrix  ♦,  then,  represents  the  sensitivities  of  the 
states  at  time  "t"  to  small  changes  in  the  initial  conditions.  It.  is 
determined  by  numerically  integrating  the  matrix  differential  equation 


*(t,tQ)  =  *(t,tQ)  =  AC t )  *(t,to),  (1-26) 

with  initial  conditions  ♦  ( tQ , tQ )  =  I,  the  6x6  Identity  matrix.  Thus, 
the  nonlinear  equations  of  motion  in  (1-11)  through  (1-13)  and  the 
matrix  equation  (1-26)  combine  to  result  in  42  first-order  differential 
equations  that  can  be  simultaneously  integrated  numerically  to 
determine  the  state  vector  and  its  associated  transition  matrix  at  any 
instant  of  time.  The  reference  trajectories  that  are  of  interest  in 
this  research  are  generated  by  a  numerical  integration  method  that  uses 
a  differential  corrections  process  developed  by  Howell  and 

[y-14] 

Pernicka.  The  orbits  Include  solar  radiation  pressure  forces  as 
formulated  by  Bell1151  specifically  for  an  orbit  associated  with  the 
interior  Lagrange  point  in  the  Sun-Earth  system.  The  numerical 
integration  routines  used  in  this  work  are  fourth-  and  fifth-order 
Runge-Kutta  formulas  available  in  the  386-Matlab  software  package. 1161 


F.  Bounded  Orbits  Near  Libration  Points 

The  computation  of  bounded  periodic  and  quasi-per iodic  orbits  in 
the  vicinity  of  libration  points  has  been  of  increasing  interest  during 
the  past  100  years.  This  section  first  discusses  the  stability  of  the 
libration  points  L.  the  CR3BP  and  the  ER3BP.  The  construction  of 
bounded  orbits  near  the  collinear  Lagrange  points  is  then  summarized. 
Finally,  the  specific  reference  trajectories  used  in  the  orbit 
determination  error  analysis  and  station-keeping  studies  in  this  work 
are  introduced. 
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1.  Stability  of  the  Libration  Points  in  the  CR3BP 


The  accomplishments  of  those  researchers  who  have  constructed 
bounded  orbits  near  collinear  libration  points  are  particularly 
significant  because  the  collinear  points  are  considered  "unstable" 
points  of  equilibrium  but  with  (only)  one  mode  producing  positive 
exponential  growth.  Bounded  motion  in  their  vicinity,  therefore,  is 
determined  by  deliberately  not  exciting  the  unstable  mode.  A  second 
mode  produces  negative  exponential  orbital  decay  and  is  also 
deliberately  not  excited.  In  the  CR3BP,  the  remaining  four  eigenvalues 
are  purely  imaginary.  The  existence  of  initial  conditions  that  result 
in  only  trigonometric  (sinusoidal)  functions  as  solutions  means  that 
the  collinear  libration  points,  while  unstable,  possess  conditional 
stability  (with  proper  choice  of  initial  conditions)  in  the  linear 

IIVI 

sense. 

The  triangular  libration  points  are  marginally  stable  in  the 

linear  sense  for  a  specific  range  of  primary  mass  ratio  in  the  CR3BP. 

Purely  imaginary  roots  in  two  conjugate  pairs  are  obtained  for  ps.0385, 

which  is  given  here  to  four  decimal  places  and  is  sometimes  referred  to 
( 18  ] 

as  Routh’ s  value.  The  mass  ratios  (listed  here  to  three  decimal 

places),  for  example,  in  the  three-body  systems  of  the  Earth-Moon 
(p  =  1.216  x  10  6),  Sun-Ear th+Moon  (p  =  3.022  x  10  6)  and  Sun-Jupiter 

-4 

(p  =  9.485  x  10  )  all  satisfy  the  mass  ratio  requirement  for  marginal 
stability  of  the  triangular  points  in  the  linearized  model.  Natural 
satellites,  such  as  the  Trojan  asteroids  or  a  moon  of  Saturn,  occupy 
linearly  stable  orbits  near  triangular  libration  points  in  their 
respective  systems. 


2.  Stability  of  the  Libration  Points  in  the  ER3BP 

Several  researchers  have  analyzed  the  stability  of  the  libration 
points  in  the  elliptic  problem,  where  both  the  mass  ratio,  p,  ?nd  the 
primary  orbit  eccentricity,  e,  influence  stability. 117-211  The 
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instability  of  the  collinear  libration  points  as  determined  in  the 
circular  problem  for  all  the  values  of  mass  parameter  persists  for  the 
elliptic  problem;  an  analysis  of  the  collinear  points  shows  instability 
for  any  combination  of  the  values  of  both  p  and  e. 

The  results  of  a  linearized  stability  analysis  regarding  the 

effects  of  eccentricity  and  mass  ratio  on  the  linear  stability  of  the 

triangular  points  have  been  published  by  Danby1201  and  then  later  by 
[21 1 

Bennett  .  Both  Danby  and  Bennett  have  numerically  generated  graphic 
depictions  of  the  linear  stability  region  in  the  p-e  plane.  For  the 
eccentricity  in  the  Sun-Earth+Moon  ER3BP,  the  value  of  ft  which  ensures 
linear  stability  is  only  slightly  less  than  Routh’s  value  (decreased  by 
approximately  one  percent).  An  interesting  aspect  of  the  p-e  stability 
region  is  that  a  range  of  values  of  p  greater  than  Routh’s  value  also 
defines  a  region  of  linear  stability  for  a  specific  range  of  values  of 
e  less  than  .3143. 


3.  Construction  of  Bounded  Collinear  Libration  Point.  Orbits 

The  initial  goal  in  the  process  of  generating  bounded  orbits  near 
a  collinear  (unstable)  libration  point  is  to  avoid  exciting  the 
unstable  mode  associated  with  the  linearized  motion.  The  meteoric  dust 
particles  that  may  be  orbiting  near  Lagrange  point  L2  in  the  Sun-Earth 
system  could  only  linger  near  that  point  if  they  arrive  with  the 
"correct"  initial  position  and  velocity  states  relative  to  L^.  The 
"correct"  initial  conditions  will  only  (primarily)  excite  the  stable 
modes  associated  with  the  linearized  motion  and  net  lor  minimally) 
excite  the  unstable  mode.  The  degree  to  which  the  unstable  mode  is 
excited  will  determine  the  length  of  time  that  the  dust  particles 
linger  near  L2_ 

The  third-order  analytic  representation  is  used  in  this  work  to 

provide  the  initial  model  for  the  trajectories.  The  method  of 

successive  approximations,  using  the  linearized  solution  as  the  first 

appi  oximation  to  the  nonlinear  orbital  path,  and  the  method  of  dual 

(6  7  22] 

time  scales  are  used  to  derive  the  third-order  result.  ’  ’  The 

method  of  successive  approximations  is  used  to  generate  an  asymptotic 
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series  in  an  appropriately  small  parameter.  (The  square  root  of  the 
eccentricity  of  the  primary  orbit,  that  is  the  orbit  of  the  Earth-Moon 
barycenter  about  the  Sun,  is  the  small  parameter  used  here. )  The 
method  of  dual  time  scales  is  used  to  convert  the  system  of  ordinary 
differential  equations  into  a  system  of  partial  differential  equations. 
In  general,  the  method  of  multiple  scales  allows  the  various  nonlinear 
resonance  phenomena  to  be  included  in  the  approximate  analytic  solution 
and  provides  a  method  to  remove  secular  terms.  (Here,  "secular"  refers 
to  terms  that  include  the  time  variable  and  is  derived  from  the  French 
"sidcle"  meaning  century.) 

The  analytic  solution  of  Richardson  and  Cary  for  the 
Sun-Ear th+Moon  ER3BP  has  been  derived  to  fourth  order,  but  the  third- 

19-14] 

order  approximation  is  found  to  be  sufficient  for  this  research. 

A  numerical  integration  algorithm,  using  a  differential  corrections 

procedure  that  is  designed  to  adjust  the  first  guess  as  obtained  from 

the  analytic  approximation,  can  then  be  used  to  numerically  generate 

the  orbit  of  interest.  A  method  developed  by  Hcwell  and  Pernicka19’141 

is  used  here  to  generate  the  orbital  paths.  Their  method  initially 

employs  the  approximate  analytic  solution  to  compute  target  points.  A 

two-level  (position  matching  then  velocity  matching),  multi-step 

differential  corrections  algorithm  is  used  to  construct  a  numerically 

integrated,  bounded  trajectory  that  is  continuous  in  position  and 

115] 

velocity.  A  solar  radiation  pressure  model  developed  by  Bell  is 
also  Incorporated  in  the  numerical  integration  procedure. 

The  method  of  Howell  and  Pernicka,  including  solar  radiation 
pressure,  uses  an  initial  analytic  guess  that  represents  a  halo  orbit 
or,  alternatively,  a  considerably  smaller  Lissajous  path.  The 
higher-order  terms  tend  to  slightly  alter  the  first-order  periodic  or 
quasi-periodic  paf‘  .  Consequently,  the  initial  target  path  for  a  halo 
orbit  will  generally  not  be  precisely  periodic.  The  two-level, 
multi-step  differential  corrections  procedure  then  adjusts  the  initial 
analytic  target  orbit  and,  therefore,  will  compute  a  halo-type  orbit 
that  is  nearly  (but  not  exactly)  periodic. 
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4.  The  Reference  Paths  Generated  for  This  Work 


Precisely  periodic  halo  orbits  exist  in  the  CR3BP.  They  also 
exist  in  the  ER3BP,  but,  in  the  ER3BP,  they  are  multiple  revolution 
trajectories  with  periods  much  longer  than  those  of  Interest  here. 
Nearly  periodic  orbits  are  more  practical  in  the  ER3BP  and  are  much 
more  likely  to  be  used  in  mission  planning;  therefore,  the  goal  here 
should  be  slightly  modified  to  be  the  comparison  of  Lissajous  and 
"halo-type"  orbits.  The  general  shapes  of  the  three-dimensional 
halo-type  and  Lissajous  orbits  can  be  seen  by  plotting  three 
orthographic  views  of  each  orbit,  using  the  tabular  data  from  the 
numerical  integration  routine.  Figure  1-4  depicts  three  orthographic 
views  of  point  plots  for  the  Lissajous  orbit  used  in  this  research. 
Figure  1-5  contains  three  orthographic  views  (on  a  slightly  different 
scale)  of  the  considerably  larger  halo-type  orbit.  (Note  that,  in 
general,  the  amplitude  ratio  for  Lissajous  trajectories  is  arbitrary. 
In  halo  orbits,  however,  constraining  the  amplitude  ratio  results  in 
equalized  frequencies  for  in-plane  and  out-of-plane  motion. )  The 
orbits  are  depicted  in  the  rotating  reference  frame  centered  at  L  . 

Both  orbits  are  clearly  not  periodic;  a  Lissajous  orbit  is  often 
called  a  quasi-periodic  path,  and  these  two  orbits  could  clearly  be 
termed  quasi-periodic  or  Lissajous  paths.  The  major  difference  between 
the  orbits  is  the  larger  size  of  the  halo-type  orbit;  however,  other 
differences  are  also  present.  The  maximum  x  and  y  excursions  of  the 
halo-type  orbit  are  approximately  four  times  as  large  as  those  of  the 
Lissajous  path.  Furthermore,  the  direction  of  motion  (clockwise  versus 
counterclockwise),  as  viewed  in  the  y-z  orthographic  depiction,  is 
different  for  the  two  orbits  used  here.  The  direction  of  motion  on  the 
halo-type  orbit  is  counterclockwise  in  the  y-z  depiction;  the  direction 
of  motion  is  clockwise  in  the  y-z  depiction  for  the  Lissajous  path. 
(Both  orbits  include  clockwise  motion  in  the  x-y  depiction. ) 

The  two  orbits  can  also  be  differentiated  in  terms  of  the 
direction  of  the  maximum  z  excursion  in  the  x-z  depiction.  If  the 
maximum  z  excursion  is  in  the  positive  z  direction,  the  orbit  can  be 
termed  a  member  of  a  "northern  family"  of  orbits.  When  the  maximum  z 
excursion  of  the  orbit  is  in  the  negative  z  direction,  the  orbit  is 
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Figure  1  -4.  Three  Orthographic  Views  of  a  Lissajous  Orbit. 
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.  indicates  the  libration  point  location. 


Figure  1-5.  Three  Orthographic  Views  of  a  Halo-Type  Orbit. 
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termed  a  member  of  a  "southern  family"  of  orbits.  In  the  x-z 
orthographic  depiction,  the  smaller  (Lissajous)  path  can  be  seen  to  be 
a  member  of  a  northern  family  of  orbits,  while  the  halo-type  orbit  is  a 
member  of  a  southern  family  of  orbits. 

Future  work  with  these  two  orbits  will  include  studies  that 
generally  require  access  to  a  nominal  path  that  is  at  least  piecewise 
smooth.  Some  method  of  curve  fitting  the  numerically  integrated  data 
must  consequently  be  investigated,  and  the  evaluation  of  various 
methods  can  add  valuable  problem  insight. 
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CHAPTER  2:  CURVE  FITTING  THE  NOMINAL  PATH 


The  reference  solution  used  in  this  research  is  generated  by 

numerical  integration  of  the  nonlinear  equations  of  motion.  In  one 

study,  an  investigation  that  used  a  consistent  dynamic  model  for  all 

comparisons,  Richardson  has  shown  that  a  slight  reduction  In  fuel 

expenditure  can  be  realized  if  a  numerically  integrated,  rather  than  an 

approximate  analytic,  nominal  path  is  computed.  The  numerical 

[9-14] 

integration  method  developed  by  Howell  and  Pernicka  is  used  here 

to  generate  a  set  of  reference  points  for  both  position  (three  states) 
and  velocity  (three  states),  relative  to  the  libratlon  point  of 
interest,  at  specified  times.  Time  series  point  plots  of  all  six  state 
variables  for  approximately  a  2-year  segment  of  a  Lissajous  orbit  are 
depicted  in  Figures  2-1  (position  states)  and  2-2  (velocity  states). 
The  method  computes  numerical  data  for  the  six  states  in  a  reference 
frame  that  is  centered  at  the  libration  point  (in  this  case  L  )  and 
that  rotates  with  the  primaries.  However,  state  estimation  techniques 
and  station-keeping  algorithms  considered  in  follow-on  research  require 
access  to  a  continuous  nominal  path,  rather  than  point  plots,  of 
acceptable  accuracy. 

The  reference  trajectory,  represented  as  a  (piecewise)  smooth 
curve,  could  be  constructed,  approximately  or  exactly,  through  the 
points  obtained  from  the  numerical  integration  routine.  The  work  here 
assumes  that  a  curve  that  passes  through  the  numerical  data  (exactly) 
is  preferred.  The  effort  required  to  generate  a  numerical  solution, 
including  forces  modeled  consistent  with  the  ER3BP  (or  even  more 
accurately  modeled  with  ephemeris  data)  would  seem  to  be  wasted  if  the 
reference  curve  deviates  too  far  from  the  numerical  data.  However,  a 
method  that  approximates  a  smooth  curve  through  the  points  is  also 
desirable;  that  is,  linear  interpolation  between  the  numerical  data 
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Figure  2-1.  Time  Series  Plots  of  Three  Lissajous  Position  States. 
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points  was  not  considered  acceptable.  In  one  study,  Per nicka1121 
that  station-keeping  costs  for  a  libration  point  orbit  were,  in  (act, 
sensitive  to  the  accuracy  of  the  curve  fit.  Clearly,  a  piecewise 
linear  curve  fit  could  not  accurately  match  the  concavity  of  the  actual 
orbital  path  between  data  points,  regardless  of  the  size  of  the  time 
steps  used  in  the  numerical  integration  routine. 

Four  methods  of  generating  a  curve  for  the  nominal  trajectory  have 
been  evaluated:  Fourier  series,  least  squares,  weighted  least  squares, 
and  cubic  splines.  The  states  associated  with  a  quasi-periodic  path 
were  thought  to  be  the  most  difficult  to  curve  fit;  therefore,  various 
Lissajous  trajectories  were  used  to  evaluate  the  curve-fitting  methods. 
One  of  these  methods  is  then  selected  to  be  used  for  the  reference 
trajectory  in  follow-on  state  estimation  and  station-keeping  research 


A.  Integrated  Fourier  Series 

Joseph  Fourier  (1768-1830)  solved  heat  flow  problems  using  the 
method  of  separation  of  variables,  and  this  application  led  to  the 
development  of  the  Fourier  series.  Leonhard  Euler  and  Alexis-Claudo 
Clairaut  (1713-1765)  had  already  expanded  some  functions  in  such  series 
representations  and  had  obtained  the  general  integral  equations  for  the 
coefficients.  Much  earlier,  the  Babylonians  had  actually  used  a 
rudimentary  form  of  a  Fourier  series  for  estimation  problems. 1361 
However,  Fourier  made  the  crucial  observation  that  every  function  could 
be  so  represented  even  if  the  function  was  not  periodic  or  not  even 

[23]  [8  22] 

continuous.  Richardson  ’  used  numerically  integrated  Fourier 

series  to  represent  the  three  position  coordinates  corresponding  to  a 
halo  orbit  (near  in  the  Sun-Earth  system).  Richardson’s  work,  and 
the  periodic  or  quasi-periodic  nature  of  the  orbits  that  require  a 
curve  fit,  made  the  initial  investigation  of  Fourier  series 
representations  seem  natural. 

The  three  position  and  three  velocity  variables  can  be  expressed 
parametrically  as 

x  =  x(t),  y  =  y(t),  z  =  z(t) ,  x  =  x(t),  y  =  y(t),  z  =  z(t),  (2-1) 
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and  the  plotted  time  history  of  the  numerical  value  of  each  variable, 

as  seen  in  Figures  2-1  and  2-2,  reveals  approximately  sinusoidal 

behavior.  (The  dot  denotes  the  derivative  with  respect  to  time. )  The 

six  states  can  be  expressed  in  separate  Fourier  series,  or  the 

expressions  for  the  three  velocity  states  can  be  computed  from  the 

first  time  derivatives  of  the  series  expressions  for  the  position 

[8  22] 

variables.  Richardson  ’  used  the  latter  method  to  find  series 

representations  for  the  velocity  components  corresponding  to  a 

libration  point  halo  orbit.  Similar  techniques  have  also  been  used  for 
a  Lissajous  path. 

Studies  by  Howell  and  Pernicka  have  used  numerically  integrated 
Fourier  coefficients  for  parameterized  series  expressions  for  the 
x=x(t),  y=y(t),  and  z=z(t)  coordinates  of  position  in  a  Lissajous 

orbit. 113,141  The  series  coefficients  are  integrated  concurrently  with 
the  equations  of  motion. 

For  example,  the  coefficients  for  the  x  position  series 
representation  would  be  computed  using 


where 


x(t)  =  — ^  +  V  [a  cos (—J—)  +  b  sin(~),  (2-2) 

n  L  1  L  1  L 


1  =  1 


and 


21 


a  = 


1 

L 


x(t)  cos(-i^)  dt,  l 


=  0,  1,2 . n  (2-3) 


2L 


b  = 


1 

L 


x ( t )  sin('~)  dt,  1=0,  1,  2, 


,  n. 


(2-4) 


The  variable  "L"  is  selected  such  that  2L  is  the  duration  of  the  orbit. 
(The  duration  of  the  orbit  is  sometimes  referred  to  in  this  work  as  the 
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"fundamental  period.")  The  variable  "n"  determines  the  number  of  terms 

retained  in  the  series;  for  instance,  when  n  =  16  and  when  a  is 

o 

assumed  to  be  zero,  the  resulting  series  representation  for  x(t)  will 
have  2n  =  32  terms.  Note  that  there  must  also  be  series 

representations  for  y ( t )  and  z(t).  (There  may  also  be  series 
representations  for  velocities  x(t),  y(t),  and  z(t),  or  time 

derivatives  of  the  Fourier  series  representations  of  the  position 
states  may  be  used,  depending  on  the  method  selected. )  The  Fourier 
series  coefficient  equations  are  numerically  inte^-ated  along  with  the 
equations  of  motion  (1-11)  through  (1-13).  When  only  the  position 
states  are  considered,  this  numerical  integration  then  includes  102 

first-order,  coupled  differential  equations  for  computation  of  the 

series,  using  32-terms  for  each  of  the  position  states. 

A  comparison  of  the  numerical  output  with  the  32-term  Fourier 
series  representations  for  the  spatial  coordinates  shows  variation  up 
to  several  thousand  kilometers.  Consistent  with  Richardson,  the  series 
representations  for  the  velocity  states  were  initially  computed  by 
taking  the  first  time  derivatives  of  the  Fourier  series  for  the 
position  coordinates.  This  method  leads  to  a  nominal  path  that  is 
different  in  position  and  velocity  from  the  results  seen  in  the 

numerical  output.  Differences  up  to  several  thousand  kilometers  in 
each  of  the  position  states  were  considered  unacceptable  and  provoked 

[14] 

further  investigation  into  curve-fitting  techniques.  As  one 

alternative,  the  number  of  terms  in  each  series  could  be  increased  and 
the  coefficients  for  the  nominal  velocity  states  could  be  numerically 
integrated  to  improve  accuracy,  but  this  technique  would  add  to  an 
already  high  computational  burden.  Also,  by  selecting  2L  to  be  the 
full  duration  of  the  orbit,  the  half-period  will  be  so  large  that  the 
first  several  coefficients  of  each  series  will  not  significantly 
contribute  to  the  series  solution.  Some  of  these  terms  could, 
consequently,  be  neglected;  however,  another  method,  such  as  least 
squares,  might  provide  closer  agreement  with  the  numerical  output  and 
could  lesson  the  computational  requirements. 
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B.  Least  Squares  Curve  Fits 


The  least  squares  approach  constructs  the  "best"  curve  through  a 
set  of  points  in  the  sense  of  minimizing  the  sum  of  the  squared  errors 
in  the  locations  of  the  points  relative  to  the  curve.  This  powerful 
method  has  had  many  contributors,  a  few  of  which  are  mentioned  here. 
The  Babylonians  set  the  stage  for  the  development  of  the  method  of 
least  squares  by  using  a  rudimentary  form  of  the  Fourier  series  for 

estimation  problems;  then  in  1632,  Galileo  Galilei  first  attempted  to 

[24] 

minimize  various  functions  of  estimation  errors.  In  1795,  Carl 

Friedrich  Gauss  invented  and  used  the  method  of  least  squares  when  he 

[23 } 

was  just  18  years  old.  Adrien-Marie  Legendre  independently 

[23’ 

invented  and  first  published  his  least  squares  development  in  1806. 

The  least  squares  approach  was  also  developed  independently  in  America 

[24] 

by  R.  Adrain  in  1808.  Gauss  employed  his  least  squares  techniques 

in  an  orbit  measurement  problem  and  developed  data  processing  methods 

125) 

for  dealing  with  random  variables.  In  order  to  help  astronomers 

locate  the  asteroid  Ceres,  Gauss  invented  a  recursive  (sequential) 
least  squares  procedure  to  handle  the  vast  calculations  that  were 

[24] 

necessary.  New  data  could  then  be  added  sequentially  rather  than 

the  problem  being  recalculated  with  the  complete  batch  of  data  when  one 
new  observation  arrived. 

The  task  of  finding  representations  for  the  numerical  output  can 
be  viewed  as  the  problem  of  fitting  six  "best"  curves  through  the  six 
given  sets  of  points:  (t,x(t)),  (t,y(t) ) ,  ( t , z ( t ) ) ,  ( t , x ( t ) ) , 

(t,y(t)),  and  ( t , z ( t ) ) .  (The  dot  denotes  the  first  derivative  with 
respect  to  characteristic  time. )  The  method  of  least  squares  can  be 
used  for  this  purpose.  Initially,  the  basis  functions  that  are  used  to 
construct  the  best  curve  through  the  data  points  must  be  selected. 
This  choice  of  basis  functions  is  usually  made  by  using  some 
hypothesized  relationship  for  the  variables  as  functions  of  time  or 
through  an  analysis  of  a  plot  of  the  points.  A  set  of  observations  of 
x(t)  .*».t  distinct  times  t  ,  t  ,  t  is  assumed.  Suppose  that  x(t) 

12  m 

can  best  be  represented  by  a  linear  combination  of  n  basis  functions. 
These  basis  functions  could  be 
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or 


1.  t, 


1,  cos(t),  sin(t),  cos(2t),  sin(2t), 


Choosing  basis  functions  that  are  orthogonal  can  lessen  the 
possible  numerical  difficulties  encountered  later  due  to  dependencies; 
however,  in  many  cases,  this  is  not  necessary.  One  orthogonal  choice 
could  be  Legendre  polynomial  (P  ,  .....  P  )  where  the  polynomial  basis 

0  n 

(1,  t,  \T .  tn)  is  orthogonal ized  using  Gram-Schmidt  methods  to 

2 

give  P  =  1,  P  =  t,  P  =  t  -  1/3 .  Another  basis  choice  might 

be  the  orthogonal  Chebyshev  polynomials  where  a  cosine  series  such  as 
(1,  cos(0),  cos(20)  =  2cos2(0)  -1,  cos(30)  =  4  cos3(0)  -  3  cos  (0),  ..) 
is  converted  using  a  substitution  "t  =  cos(0)"  to  give  T  =  1,  T  =  t, 
T  =  2  t2  -1,  T  =  4  t3  -  3  t .  A  third  choice  of  ort  ogonal 

2  3 

basis  functions  could  be  the  sinusoids  1,  cos(t),  sin(t),  cos(2t), 
sin(2t) . 

In  general,  the  method  can  be  illustrated  by  assuming  a  set  of  n 
basis  functions  denoted  by  1,  V(t),  V(t),  V(t),  .  V  (t). 

12  3  n-1 

Then  an  expression  for  x(t)  can  be  computed  using  the  basis  functions 
V4 ( t )  and  the  scalar  coefficients — as  yet,  unknown — denoted  as  b$  for 
l  =  1,  2,  3 . n-1; 

x(t)  =  b  +  b  V  (t)  +  .  +  b  V  (t)  -  error. 

0  11  n-1  n-1 


The  term  "error"  is  included  because  the  approximation  is  not,  at  the 
outset,  assumed  to  be  exact.  With  m  observations  and  with  the  error 
denoted  by  "e",  we  can  write  m  separate  equations: 
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b  +  b  V  (t  )  + 
o  111 


+  b  V  (t  )  -  e  =  x(t  )  =  x 

n-l  n-1  11  11 


b  +  b  V  (t  )  +  .  +  b  V  (t)-e  =  x(t  )  =  x 

0  112  n-l  n-l  2  2  2  2 


b  +  b  V  (t  )  +  .  +  b  V  (t  )  -  e  =  x(t  )  =  x  (2-5) 

0  113  n-l  n-l  3  3  1  3 


b  +  b  V  (t  )  + .  +  b  V  (t)-e  =  x(t  )  =  x 

0  11m  n-l  n-l  mm  mm 


This  system  (2-5)  of  m  equations  in  n  unknowns  can  be  represented  in  a 
matrix  equation  as 


A  b  -  e  =  x 


(2-6) 


where 


1  V  ( t  )  V  (t  )  .  V  (t  ) 

11  21  n-l  1 


1  V  (t  )  V  (t  )  .  V  (t  ) 

1  2  2  2  n-l  2 


A  =  1  V  (t  )  V  (t  )  .  V  (t  )  ,  (2-7) 

1  3  2  3  n-l  3 


1  V  (t  )  V  (t  )  .  V  (t  ) 

1  m  2  m  n-l  m 
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In  general,  least  squares  problems  are  overdetermined  (with  m  >  n);  an 
exact  solution  to  the  overdetermined  system  A  b  =  x  is  possible  only  if 
the  error  vector,  e,  is  the  zero  vector  (all  m  equations  hold  exactly). 
Generally,  an  exact  solution  is  not  possible.  There  are  many  methods 

r  26] 

that  can  be  used  to  solve  an  equation  such  as  (2-6);  Skelton 
discusses  several  of  the  methods. 

The  least  squares  method  minimizes  the  sum  of  the  squared  error 
terms.  Since  any  particular  error  term  could  be  positive  or  negative, 
the  sum  of  the  errors  could  be  very  small  while  individual  errors  were 
still  quite  large.  Consequently,  the  sum  of  the  squared  errors  is  the 
cost  function  that  is  minimized.  We  minimize  P  where 


P  =  £  e*  =  eTe.  (2-11) 

1  =  1 


The  overbar  denotes  a  vector,  and  the  "t"  indicates  a  transpose.  Using 
equation  (2-6)  in  equation  (2-11),  the  cost  function  simplifies  to 


P  =  bTATA  b  -  2  xTA  b  +  xTx. 


(2-12) 


Minimizing  equation  (2-12)  with  respect  to  b  results  by  taking  the 
derivative  of  the  scalar  P  with  respect  to  the  vector  b: 

—  =  ATAb  +  (bTATA)T  -  2  (xTA)T  =  0.  (2-13) 

3b 


Note  that  the  term  xTx  in  equation  (2-12)  is  assumed  to  be  a  constant 
and  does  not  appear  in  the  derivative  (2-13).  Simple  matrix  algebra 
then  leads  to  the  solution 


b  =  (ATA)-1ATx. 


(2-14) 
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where  b  is  used  to  denote  the  least  squares  estimate  of  the  "best" 
vector  of  basis  coefficients. 

Equation  (2-14)  is  a  depiction  of  the  solution;  however,  there  are 
numerous  solution  methods  that  are  more  numerically  stable  than 
inversion  of  such  a  large  system  as  A  A  in  (2-14).  For  instance,  when 
the  columns  of  A  are  not  independent  (as  can  occur  when  the  basis 
functions  are  not  orthogonal),  A  A  will  not  be  invertible,  and  the 

T 

psuedoinverse  of  A  A  can  be  used  to  complete  the  least  squares 

1271  T  T  -1 

approach.  Even  when  A  A  is  invertible,  the  notation  (A  A)  is 

included  only  to  depict  a  solution  for  the  unknown  coefficients. 

Generally,  a  more  computationally  efficient  and  numerically  stable 

method  than  matrix  inversion  is  used.  A  method  such  as  QR 

factorization  using  Gram-Schmidt  methods  or  the  use  of  a  set  of 

T  [27] 

Householder  transformations  is  more  efficient  than  inverting  A  A. 

In  order  to  briefly  illustrate  the  QR  decomposition  method,  the  matrix 
T 

A  A  is  written  as 


ATA  =  QR 


—  1  T 

where  Q  is  an  orthogonal  matrix  such  that  Q  =  Q  and  and  R  is 
triangular  and  easily  inverted.  The  least  squares  solution  derived 
from  equation  (2-13)  is 


ATA  b  =  ATx 


Using  the  QR  decomposition,  this  equation  can  then  be  written  as 


(QR)b  =  ATx. 
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The  solution 


A  -IT  T— 

b  =  R  QAx 

involves  inversion  of  only  a  diagonal  matrix  R. 

For  this  work,  sinusoidal  basis  functions  were  the  logical  choice 
for  the  states  (see  Figures  2-1  and  2-2).  Using  the  x  coordinate  as  an 
example,  one  row  in  the  least  squares  formulation  for  the  time  step  m 
can  be  written  as 


a  +a  cos 
o  1 


filn  Cm 
+a  COSI  — 
n  '  L 


(2-15) 


where  "L"  is  the  half-period  selected  for  the  computation.  The  number 
of  basis  functions  (2n+l)  can  be  chosen  to  pr >duce  the  most  accurate 
representation,  as  well  as  limiting  the  computational  burden.  If  101 
basis  functions  we-e  desired,  the  "last"  basis  function  would  be 
sinf — ,  and  its  (yet  to  be  determined)  coefficient  would  be  b  . 
The  rows  representing  all  integration  time  steps  can  be  combined  in 
matrix  form  to  produce  the  matrix  equation 


A  b  =  x 


(2-16) 


where  A  is  a  matrix  with  m  rows,  one  for  each  of  the  m  time  steps,  and 
with  n  columns  corresponding  to  the  number  of  trigonometric  basis 
functions  selected.  The  vector  b  contains  the  coefficients  for  the 
basis  functions.  The  column  vector  x  has  m  entries  for  the  x 
coordinate,  one  for  each  of  the  m  time  steps.  (Note  that  if  curve  fits 
for  all  six  states  were  necessary,  then  this  method  would  be  applied 
six  times  with  the  right  hand  side  of  equation  (2-16)  being  x  then  y 

A 

and  so  forth.  The  coefficients  in  b  determined  for  each  state  variable 
would  consequently  be  different. ) 

A 

The  least  squares  approach  solves  for  the  "best"  estimate,  b,  of 

[28] 

the  vector  b  by  minimizing  the  sum  of  the  squared  errors.  The 
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estimate  b  is  calculated  such  that  x  =  Ab  differs  from  x  by  some  error 
vector  e,  where 


e  =  x  -  x, 


(2-17) 


and  the  least  squares  estimate  that  results  is 


b  =  (AtA)'Vx. 


(2-18) 


The  scalar  entries  in  the  vector  b  are  then  the  Fourier  coefficients 
for  the  basis  functions  selected  for  use  in  the  formulation.  This 

approach  to  calculating  the  coefficients  of  the  sinusoidal  basis 
functions  has  the  advantages  of  computational  efficiency  and  great 
flexibility  in  the  choices  of  fundamental  period,  2L,  and  the  number  of 
basis  functions,  2n+l,  retained  in  the  series. lZ9]  By  choosing  basis 
functions  that  are  orthogonal,  possible  numerical  difficulties  due  to 
dependencies  can  be  eliminated,  and  the  independent  contribution  of 
each  basis  function  can  be  judged  by  the  relative  size  of  its 
coefficient. 

Two  criteria  are  used  here  to  compare  the  various  options  for  the 
fundamental  period  and  the  number  of  basis  functions.  The  value  of  the 
sum  of  the  squared  errors  (eTe)  can  be  calculated  in  order  to  compare 
various  choices  of  n  and  L.  The  element  of  the  error  vector  e  with  the 
largest  magnitude  provides  an  indicator  of  the  maximum  error  for  a 
specific  choice  of  n  and  L,  and  it  is  denoted  as  e  .  The  scalars 

max 

e  can  also  be  used  as  a  criteria  for  evaluating  various  combinations 

max 

of  n  and  L.  The  error  vector  is  calculated  by  noting  that 


e  =  x-  x  =  Ab-x  =  A(ATA)"Vx  -  x  =  (A(ATA)'V-  I]  x  (2-19) 


The  identity  matrix  is  denoted  as  I.  The  value  e  is  then  defined  as 

max 

the  entry  in  the  vector  e  that  has  the  largest  absolute  value.  In  one 
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evaluation  of  this  method,  a  Llssajous  orbit  of  approximately  5  years 
duration  was  used. 

The  quasi-periodic  nature  of  the  Ll.'sajous  orbit  allowed  some 

flexibility  in  the  choice  of  "period"  for  modeling.  In  one  trial,  a 

half-period  of  L  =  87.5  days  was  roughly  estimated  from  an  analysis  of 

the  plot  of  x(t)  versus  time.  Other  values  of  L  close  to  87.5  days 

were  also  used  in  separate  simulations,  and  the  criteria  involving 

determination  of  eTe  and  e  were  evaluated.  A  half-period  equal  to 

max 

the  duration  of  the  orbit  (L  =  1741.5  days)  was  also  used  and  was 
expected  to  minimize  the  error  vector  magnitude  for  "large"  values  of 

[29] 

n.  Over  numerous  trials,  the  best  curve  fit  was  obtained  for  the 
combination  L  =  1741.5  days  and  n  =  101.  The  results  of  several 
simulations  are  presented  in  Table  2-1. 


TABLE  2-1:  COMPARISON  OF  LEAST  SQUARES  CURVE  FITS 


n 


L  =  87.5  Days 


L  =  1741.56386  days 


e 

max 


(km) 


eTe  (km2) 


e  (km)  eTe(km2) 

max 


101 

6.4866 

X 

)0J 

4.9455 

X 

10 

2.7262 

X 

io2 

1.5307 

X 

10 

81 

6.4805 

X 

103 

4.9455 

X 

109 

2.7095 

X 

102 

1 . 7092 

X 

10 

61 

6.4753 

X 

io3 

4.9456 

X 

109 

1 . 3447 

X 

103 

6.6568 

X 

10 

41 

6.4820 

X 

103 

4.9457 

X 

109 

3.0322 

X 

103 

1.6871 

X 

io' 

31 

6.4750 

X 

103 

4.9457 

X 

109 

6.8196 

X 

103 

3.8327 

X 

io' 

21 

6.4930 

X 

103 

4.9459 

X 

109 

6.9039 

X 

103 

1 . 0044 

X 

10 

11 

6.5635 

X 

103 

4.9495 

X 

109 

6.3270 

X 

104 

1 . 2776 

X 

10 

3 

7.5111 

X 

103 

5.4017 

X 

109 

5.9419 

X 

io4 

1 . 3569 

X 

10 

Clearly,  using  101  basis  functions  for  x(t),  y(t),  z(t),  x(t), 
y(t),  and  z(t)  could  be  computationally  expensive.  Implementing  this 
method  with  n  =  101  still  resulted  in  maximum  errors  (e  )  in  the 

max 

range  of  hundreds  of  kilometers  in  the  x  estimate.  Errors  of  such  a 
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magnitude  provided  improved  accuracy  over  a  32-term  series  with 
numerically  integrated  Fourier  coefficients,  yet  it  was  desired  to 
further  reduce  the  errors  and  the  computational  burden,  if  possible. 
Therefore,  the  third  curve-fitting  method  attempted  was  weighted  least 
squares. 


C.  Weighted  Least  Squares 

The  method  of  weighted  least  squares  allows  the  introduction  of 
weighting  factors  for  the  errors.  The  "most  important"  errors  are 
weighted  more  heavily  (relative  to  the  other  errors),  and  the  weighted 
sum  of  the  squared  errors  is  minimized.  Gauss  developed  and  first  used 

[27  30] 

weighted  least  squares.  ’  He  selected  a  weighting  scheme  for  the 
errors  that  allowed  his  weighted  least  squares  procedure  to  provide  a 
linear,  unbiased,  minimum  variance  estimate,  and  these  properties  wiil 
be  used  in  the  derivation  of  a  state  estimation  algorithm  described 
later. 

Regardless  of  how  it  is  determined,  the  weighting  matrix,  W,  is 

multiplied  by  the  error  vector,  e,  and  the  quantity  [We]T[We]  is 
[28] 

minimized.  This  minimization  leads  to  the  estimate 


c  ^  (aVwa)-Wwx. 


(2-20) 


Several  weighting  schemes  were  considered.  Comparison  of  the  least 
squares  curve  fits  with  the  numerical  output  revealed  that  the  poorest 
agreement  of  the  curve  with  the  numerical  data  occurred  in  the  ranges 
where  the  x(t)  position  coordinate  was  near  its  maximum  value.  In  one 
trial,  therefore,  the  errors  associated  with  data  points  for  which  the 
values  of  x(t)  were  near  their  maximum  magnitude  were  given  a 
relatively  large  weight.  The  goal  was  to  force  the  least  squares  curve 
through  these  points.  A  second  trial  involved  a  two-step  process. 
Initially,  an  unweighted  least  squares  procedure  was  completed.  The 
elements  of  the  error  vector  (and  a  function  of  e  for  some  trials)  were 
used  in  step  two  as  he  diagonal  entries  of  the  weighting  matrix.  The 
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sum  of  the  squared  errors  (eTe)  and  the  value  of  e  did  not 

max 

significantly  decrease  by  using  these  various  weighted  least  squares 
approaches.  In  the  continuing  search  to  improve  the  curve  fit,  the 
fourth  method  used  was  based  on  a  piecewise  fit  that  forced  agreement 
at  the  data  points. 


D.  Piecewise  Curve  Fits 

Early  work  using  a  least  squares  method  clearly  indicated  that  an 
approximate  (that  is,  the  "best"  approximate)  curve  fit  is  computed 
based  on  specific  criteria.  When  segments  of  about  90  days  were  used 
for  least  squares  curve  fits,  excellent  results  were  obtained;  this 
directly  led  to  the  idea  of  using  piecewise  curve  fit  methods  for  the 
data  spanning  a  complete  trajectory.  Curve  fitting  using  local  cubic 
polynomials  computed  in  a  "spline"  algorithm  was  the  final  method 
investigated. 

Although  tne  antecedents  of  splines  appeared  in  the  technical 
literature  more  than  4  decades  ago,  it  is  generally  agreed  that  the 
formal  birth  of  splines  was  in  the  1940s,  and  the  development  is 
closely  associated  with  Professor  I.  J.  Schoenberg.  The  simplest 
type  of  spline  is  a  set  of  straight  lines  connecting  consecutive  data 
points.  Other  types  of  splines  might  be  nonlinear  with  rational 
functions,  exponentials,  algebraic  functions,  or  even  some  piecewise 
combination  of  these  functions  used  for  a  specific  problem.  The  cubic 
spline  approximates  the  position  of  a  flexible  thin  beam,  or 
draftsman’s  spline,  passing  through  all  the  points  of  a  given  data  set. 
In  this  sense,  it  has  been  found  that  a  particular  type  of  cubic 
polynomial  spline  nearly  minimizes  the  strain  energy  over  all  function 

[32J 

curves  passing  through  the  given  data  set.  This  was  the  reason 

[32  33] 

behind  Schoenberg’s  choice  of  the  word  "spline."  ’ 

A  cubic  polynomial  between  any  two  data  points  has  four  degrees  of 
freedom  (four  coefficients).  So,  for  a  cubic  spline  with  m  segments, 
the  total  possible  degrees  of  freedom  would  be  4m.  Boundary  conditions 
for  the  spline  segments  can  greatly  simplify  the  solution.  Position 
matching  at  all  m  data  points,  slope  continuity,  curvature  continuity, 


42 


and  end  point  slope  or  curvature  boundary  conditions  result  in  a  system 
of  equations  with  a  tri-diagonal  (m-2)x(m-2)  coefficient  matrix.  The 
resulting  matrix  equation  can  be  speedily  solved  and  stored  as  an  array 
of  break  points  (intersections)  and  an  array  of  associated  polynomial 
coefficients,  both  available  for  easy  access.1321 

A  simple  example  can  illustrate  both  the  computation  of  the  spline 
coefficients  and,  in  particular,  the  use  of  boundary  conditions. 
Suppose  a  set  of  5  points  is  given: 

(tj.xftj))  =  (t^Xj) 

(t2,x(t2))  =  (t2>x2) 

(t3,x(t3))  =  (t3,x3) 

(t4>x(t4))  «  (t4,x4) 

( t5 ,  x(  t5 )  )  =  (t^.Xj.) 

It  is  assumed  that  these  points  are  both  the  given  set  of  five  data 
points  and  the  points  of  intersection  of  four  cubic  polynomials.  This 
assumption  is  not  in  general  necessary  and  is  made  here  only  for 
simplicity.  Each  segment  has  a  polynomial  of  the  form 

x^t)  =  ai(t-ti)3+  bi(t-t()2+  cj(t-ti)+  d^  for  1=  1,  2,  3,  4.  (2-21) 

The  first  and  second  derivatives  with  respect  to  time  are  then  given  by 

x ( t )  =  3aj ( t— 1 1 )2  +  2  b  (t-t  )  +  c  (2-22) 

and 
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(2-23) 


x  ( t )  =  6at  ( t— t ^ )  +  2 


For  notational  simplicity,  if  (tj+J  -  t^)  =  hj  and  xCt^  =  then  for 
segment  1  with  endpoints  ( t ^ ,  x^  and  (t^,  x^) 


X|  =ai(Vt|)3  +  b,(t|-t|)2  +  ci(t«-t|)  +d. 


(2-24) 


x  =  3a  (t  -t  )2  +  2b  (t  -t  )  +  c 

l  ill  ill  l 


(2-25) 


x  =  6a  (t  -t  )  +  2b 
l  ill  l 


=  2b  =  S  , 
l  i 


(2-26) 


x  =  a  (t  -t  )3  +  b  (t  -t  )2  +  c  (t  -t  )  +  d 

l+i  l  l+i  i  l  t+i  l  l  l+i  l  l 


=ah3  +  b  h2  +  ch  +  d, 
ii  ii  ii  i 


(2-27) 


x  =  3a  (t  -t  )2  +  2b 
l+i  11+11  l 


t  )  +  c 

i  l 


=  3a  h2  +  2b  h  +  c  , 

ii  ii  i 


(2-28) 


and 
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(2-29) 


x  =  6a  (t  -t  )  +  2  b 
l+i  li+il  l 

=  6a  h  +  2b  =  S 

ll  ll+i 

Using  equation  (2-26), 

b  =  S  /2\  (2-30) 

l  i 

using  equations  (2-26)  and  (2-29), 

=  (S  -  S  )/6h  ;  (2-31) 

finally,  using  equation  (2-24), 

(2-32) 

Then,  placing  these  identities  into  equation  (2-27),  the  solution 

c  =  (x  -x  )/h  -  (2h  S  +h  S  )/6  (2-33) 

l  1+111  i  l  l  l+i 

is  derived.  Using  the  condition  that  the  slopes  from  each  side  at  all 
the  interior  points  (2,  3,  and  4)  must  match  and  modifying  equation 
(2-28)  for  the  i-i  point, 
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(2-34) 


=  3a  (h  )  +  2b  h  +  c 
l-i  l-i  l-i  l-i  1-1 


Expressions  for  a(  j,  and  can  be  readily  obtained  from 

equations  (2-31),  (2-30)  and  (2-33),  respectively.  These  expressions 
can  then  be  substituted  into  equation  (2-34),  and  the  resulting 
expression  can  be  equated  to  (2-33).  An  equation  in  the  unknowns  , 
S  ,  S  and  in  the  given  values  x  ,  x  ,  h  ,  and  h  is  then 

i  l+i  6  l’  l-i  i-i*  l 

obtained  (after  some  algebra): 
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Boundary  conditions  are  used  to  find  and  Sg  for  this  5-point 
problem.  These  conditions  could  be  =  Sg  =  0;  this  choice  results  in 
what  is  called  a  "natural"  spline.  Alternatively,  the  "not-a-knot"  end 
conditions  are  present  if  the  end  point  derivatives  are  not  specified 
and  the  required  polynomial  pieces  are  reduced  from  n-1  to  n-3.  For 
the  following  example,  the  end  point  second  derivatives  are  chosen  such 
that  S  =  S  and  S  =  S  .  The  final  result  is  an  (n-2)  x  (n-2)  or,  in 

12  s  4 

this  case,  a  3  x  3  tri-diagonal  system: 

3h+2h  h  0  1  f  S  I  xl 

1  2  2  2  2 

h  2h  +2h  h  S  =  x  .  (2-36) 

2  2  3  3  3  3 

0  h  2h  +3h  S  x 

3  3  4  J  L  4  J  4J 

Equation  (2-36)  can  be  written  in  more  compact  form  as 

A  S  =  x 

Solving  for  the  unknowns  S^,  S^,  and  then  allows  solution  of 
equations  (2-31)  through  (2-33)  for  the  coefficients  a^,  b^,  c^,  and 
for  each  of  the  polynomial  segments.  The  cubic  spline,  with  curvature 
continuity,  can  then  be  stored  as  an  array  of  break  points 
(intersections)  and  an  array  of  associated  polynomial  coefficients  for 
easy  access  and  interpolation. 

The  boundary  conditions  at  points  1  and  n  could  have  been  altered 
to  stipulate  the  s.opes  at  these  endpoints.  When  these  slopes  are 
known  due  to  some  prior  knowledge  of  the  data  set,  the  resulting  cubic 
spline  with  curvature  continuity  and  with  slope  boundary  conditions  is 
called  a  "complete"  cubic  spline.  It  has  been  determined  that  a 
complete  cubic  spline  interpolant  approximately  minimizes  the  strain 
energy  over  all  function  curves  passing  through  the  given  data  set. 
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Following  a  review  of  cubic  spline  interpolants  and  the  use  of 

some  of  these  splines  to  fit  a  representative  data  set,  a  method 

[’’41 

devised  by  Hiroshi  Akima  deserves  special  mention.  The  Akima 
spline  minimizes  extraneous  points  of  inflection  and  provides  a  smooth 
curve  that  closely  resembles  one  that  could  be  drawn  manually  through 
the  data.  The  Akima  spline  interpolation  routine  produced  excellent 
results  for  the  state  coordinates  associated  with  a  Lissajous 
trajectory.  It  was  exact  through  all  modeled  points,  consistent  with 
its  formulation. 

The  most  unique  characteristic  of  the  Akima  spline  is  that  it 

combats  "wiggles"  in  the  interpolant  curve  by  minimizing  extraneous 

points  of  inflection.  In  this  way,  the  Akima  spline  provides  a  smooth 

curve  that  closely  resembles  one  that  would  be  drawn  manually  through 

the  data  set.  Several  mathematical  methods  for  interpolating  a  single 

valued  function  result  in  a  curve  that  is  very  different  than  the 

smooth  curve  that  would  be  drawn  manually.  Methods  that  place 

polynomial  or  another  type  of  curve  through  a  data  set,  and  even  some 

[34] 

splines,  generally  produce  unnatural  points  of  inflection.  It  can 

be  shown  that  there  is  a  unique  polynomial  of  degree  n  through  n+1 
points  or  a  unique  trigonometric  polynomial  with  n  terms  through  n-1 

[351 

points.  Figures  2-3,  2-4,  and  2-5  illustrate  the  use  of  a  ten 

degree  polynomial,  11-term  sinusoidal  curves,  and  spline 
approximations,  respectively,  passing  through  the  11  points  listed  in 

[34) 

Table  2-2.  Extra  "wiggles"  (points  of  inflection)  are  obvious  in 

Figures  2-3  and  2-4.  The  spline  approximations  in  Figure  2-5  are  not 
Akima  splines,  and  they  do  exhibit  an  additional  point  of  inflection 
near  t=6. 


TABLE  2-2.  LIST  OF  POINTS  USED  IN  THE  COMPARISONS  OF  METHODS 


t 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

x(t) 

10 

10 

10 

10 

10 

10.5 

15 

60 

85 
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The  Dependent  Variable  x(t) 


Using  a  10-Degree  Polynomial  Fit 


Figure  2-3.  Using  a  10-Degree  Polynomial  Curve  Fit. 
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The  Dependent  Variable  x(t)  The  Dependent  Variable  x(t) 


Using  an  1  i -Term  Sine  Series  Fit 


Using  an  1 1-Term  Cosine  Series  Fit 


Figure  2-4.  Using  a  Sinusoidal  Fit. 
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The  Dependent  Variable  x(t)  The  Dependent  Variable  x(t) 


Using  a  Cubic  Spline  Curve  Fit  with  "Not-A-Knot"  End  Conditions 


Figure  2-5.  Using  a  Cubic  Spline  Curve  Fit. 
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The  Akima  spline  relies  on  locally  determined  slopes  using 

geometrical  relationships  at  each  junction  point.  A  spline  in  general 

removes  the  global  dependence  of  a  polynomial  on  some  remote  local 

properties:  if  the  function  to  be  approximated  by  a  polynomial  is 

badly  behaved  anywhere  in  the  interval  of  approximation,  then  the 

approximation  is  poor  everywhere.  The  Akima  spline  goes  one  step 

further,  beyond  removing  global  dependence,  by  also  calculating  local 

slopes  from  the  data  set  and  matching  these  at  each  junction  point. 

The  method  uses  five  points  P  ,  P  ,  P  ,  P  ,  P  selected  from 

K  1-2’  l-i ’  1’  i+i  1+2,  , 

1 34] 

the  data  set  to  calculate  the  slope  at  point  P  .  Boundary 

conditions  are  provided  by  estimating  two  more  points  on  each  end  of 


the  data  set  and  then  calculating  the  slopes  at  the  end  points.  For 


every  five  points, 

P  F.  P"P  ,  and 
1-1  1  1  1+1 


the  slopes  of  the  four  line 

P  P  are  denoted  by  M  , 
1+1  1+2  1 


segments 


P  P  , 
1-2  1-1 

and  M. , 

4 


respectively.  The  geometrical  relationships  used  by  Akima  lead  to  an 

[  34] 

expression  of  the  slope  at  point  P  (denoted  as  t ^ : 


When  the  denominator  of  t  is  zero,  this  slope  at  point  P^  is  given  by 


tj  =  (1/2) (M2  +  M3) 


Thus,  position  and  slope  are  specified  at  each  data  point,  and  between 
any  two  data  points,  four  conditions  are  present.  A  third-order 
polynomial  can  then  be  uniquely  determined. 

The  Akima  algorithm  and  several  other  spline  interpolation 
routines  produced  excellent  results  for  the  coordinates  of  a  libration 
point  Lissajous  trajectory  and  required  very  little  computational 
memory.  While  the  Akima  spline  in  general  produces  a  curve  fit  with  a 
minimal  number  of  points  of  inflection,  the  reference  trajectories  used 
in  this  effort  are  relatively  easy  to  curve  fit  with  a  either  a 
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"natural"  or  "not-a-knot"  cubic  spline.  Routines  available  in 
132] 

386-Matlab  are  used  primarily  in  this  research  effort. 

Several  differing  spline  routines  generated  excellent  curve  fits. 
The  computation  incorporates  forced  agreement  at  all  modeled  data 
points;  however,  analysis  of  the  curve  fit  between  data  points  is  then 
necessary.  That  is,  it  would  be  interesting  to  question  how  accurately 
the  cubic  spline  models  the  actual  orbit  between  data  points.  A 
typical  numerical  data  set  for  a  800-day  Lissajous  trajectory  was 
selected  to  evaluate  the  shape-preserving  nature  of  the  spline.  Every 
other  point  of  the  data  set  was  removed  and  stored  for  later 
comparison.  An  Akima  spline  was  then  computed  through  the  remaining 
points.  This  spline  was  used  to  predict  the  locations  of  the  points 
that  had  been  removed,  and  the  spline  prediction  was  compared  to  the 
actual  values.  The  results  showed  errors,  when  using  an  Akima  cubic 
spline,  of  less  than  one  kilometer  for  the  interior  spline  points. 
Similar  tests  using  other  cubic  spline  routines  also  showed  errors  of 
consistently  less  than  one  kilometer;  the  results  are  summarized  in 
Table  2-3  and  Figure  2-6  .,  In  most  cases,  the  difference  between  the 
spline-computed  positions  and  the  positions  cf  the  removed  points  were 
less  than  .5  kilometers.  The  difference  in  accuracy  delivered  by  the 
Akima  or  a  more  general  cubic  spline  is  much  better  (in  meters)  if  the 
corresponding  numerical  data  is  available  for  time  steps  of 
approximately  2  days  or  less. 


Table  2-3.  Curve  Fit  Accuracy  Using  Cubic  Splines 
for  the  x(t)  Coordinate  of  a  Lissajous  Path. 


Orbit  Duration 

Average  integration  time  step 
Average  change  in  x(t)  during  the  time  step 
Minimum  x(t)  coordinate  value 
Maximum  x(t)  coordinate  value 


802. 

2. 

4255. 

388. 

56078. 


Maximum  error  at  unmodeled  points  using  natural  spline 
Maximum  error  at  unmodeled  points  using  not-a-knot  spline 
Maximum  error  at  modeled  points  using  natural  spline 
Maximum  error  at  modeled  points  using  not-a-knot  spline 


9063  days 
8400  days 
0000  km 
2521  km 
0000  km 
6182  km 
6181  km 
0  km 
0  km 
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Using  a  Cubic  Spline  with  "Natural"  End  Conditions 


Using  a  Cubic  Spline  with  "Natural"  End  Conditions 


Figure  2-6.  Histograms  of  Errors  Caused  by  Curve  Fitting 
a  Lissajous  Path  Using  Cubic  Splines  when 
Every  Other  Point  is  Unmodeled. 
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The  cubic  spline  has  proven  to  be  a  useful  and  accurate 
interpolating  polynomial  for  the  coordinates  of  a  L.issajous  tra jeclot v 
about  the  interior  libration  point.  Fourier  series  and  least 
squares  approaches  using  trigonometric  polynomials  could  not  match  its 
accuracy;  computation  and  storage  do  not  seem  to  present  a  problem 
when  cubic  splines  are  used  to  represent  the  nominal  trajectory. 

Cubic  splines,  including  the  Akima  spline,  were  used  in  the 
subsequent  work  for  the  computation  of  the  numerical  entries  in  the 
Jacobian  matrix  as  it  was  evaluated  along  the  nominal  path,  and  for  the 
calculation  of  residuals  relative  to  the  reference  trajectory.  When 
range  and  range-rate  measurements  are  simulated,  the  nonlinear 
measurement  equations  are  commonly  linearized  about  the  nominal  path. 
The  cubic  spline  can  then  also  be  used  to  compute  the  time-varying 
entries  in  the  resulting  measurement  matrix.  The  Akima  cubic  spline 
software  is  available  in  the  International  Mathematical  and  Statistical 
Libraries  (IMSL)  by  using  the  routines  CSAKM  (to  compute  the  spline) 
and  CSVAL  (to  evaluate  the  spline  at  a  specified  abscissa  value).  A 
useful  cubic  spline  program  is  also  available  on  386-Matlab  by 
using  the  SPLINE  (computation)  and  PPVAL  (evaluation)  routines.  The 
routines  CSAPN  and  CSAP1  can  also  be  used  in  Matlab  to  construct 
natural  and  not-a-knot  cubic  splines,  respectively.  Spline 

representations  can  now  be  incorporated  in  both  the  error  analysis  and 
station-keeping  studies  in  future  research. 


J 
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CONCLUSION 


The  restricted  three-body  problem  is  an  important  and  interesting 
area  of  research.  The  Sun-Ear th+Moon  ER3BP,  in  particular,  is 
currently  an  area  of  vital  research  attention.  The  use  of  highly 
accurate,  numerically  integrated  nominal  trajectories,  coupled  with  the 
need  of  a  continuous  representation  of  the  path,  necessitates  the 
investigation  of  curve  fit  methods..  Other  researchers  have  used 
approximate  Fourier  series  representations  to  curve  fit  the  numerically 
integrated  libration  point  orbits.  The  use  of  such  an  approximate 
curve  near  an  accurately  determined,  numerically  integrated  path  seems, 
in  part,  to  be  wasted  numerical  integration  effort.  Other  researchers 
working  on  European  Space  Agency  missions,  have  used  tabular  look-up 
methods  (linear  interpolation)  as  a  curve  fit  method.  This 
approximation  ignores  the  curvature  of  the  data,  regardless  of  the 
duration  of  the  time  steps  used  in  the  numerical  integration.  The 
curve  fitting  using  cubic  splines  exactly  models  every  data  point  from 
the  numerical  integration  routine,  accurately  models  the  curvature 
between  the  data  points,  and  is  computationally  efficient.  Simulations 
used  to  compute  the  state  transition  matrix  and  linearized  measurement 
equations  and  to  calculate  expected  tracking  errors  and  the  cost  of 
maintaining  the  spacecraft  near  the  unstable  orbits  can  now  also 
incorporate  cubic  spline  curve  fits.  This  ,  in  turn,  will  improve  the 
accuracy  of  further  studies  conducted  in  relation  to  future  spacecraft 
missions. 
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